import numpy as np
from scipy import stats

def find_min_n_for_p(p0=0.1, alpha=0.05, delta_start=0.001, delta_end=0.05, delta_step=0.001, max_n=10000):
    z_alpha_2 = stats.norm.ppf(1 - alpha/2)
    print(f"z_(α/2) = {z_alpha_2:.4f}")
    print(f"{'p':>8} {'min_n':>8} {'示例k':>8} {'p̂':>8} {'SE':>10} {'z':>8}")

    delta = delta_start
    while delta <= delta_end:
        p = p0 + delta
        n = 10
        found = False
        while n <= max_n:
            # # k ~ B(n, p)
            # k = np.random.binomial(n, p)
            # # 样本次品率
            # p_hat = k / n
            k = n * p
            p_hat = p
            # 标准误差
            SE = np.sqrt(p_hat * (1 - p_hat) / n) if n > 0 else 0
            # z统计量
            z = (p_hat - p0) / SE if SE > 0 else 0
            # 拒绝条件
            if p_hat - z_alpha_2 * SE > p0:
                print(f"{p:8.4f} {n:8d} {k:8.0f} {p_hat:8.4f} {SE:10.6f} {z:8.4f}")
                found = True
                break
            n += 1
        if not found:
            print(f"{p:8.4f} {'未找到':>8}")
        delta += delta_step

def calculate_error_range(p0=0.1, alpha=0.05, p_values=[0.107, 0.11, 0.12, 0.13, 0.14, 0.15], trials=1000):
    """
    计算确定性方法和随机方法之间的误差
    """
    z_alpha_2 = stats.norm.ppf(1 - alpha/2)
    print(f"=== 误差分析 (α = {alpha}, p₀ = {p0}) ===")
    print(f"{'p':>8} {'理论n':>8} {'随机n_avg':>10} {'随机n_std':>10} {'误差%':>8} {'误差范围':>12}")
    
    for p in p_values:
        # 理论n值（确定性方法）
        n_theoretical = 10
        while n_theoretical <= 10000:
            p_hat = p
            SE = np.sqrt(p_hat * (1 - p_hat) / n_theoretical)
            z_theoretical = (p_hat - p0) / SE if SE > 0 else 0
            if p_hat - z_alpha_2 * SE > p0:
                break
            n_theoretical += 1
        
        if n_theoretical > 10000:
            print(f"{p:8.3f} {'未找到':>8}")
            continue
        
        # 随机方法模拟
        n_random_list = []
        for _ in range(trials):
            n_random = 10
            found_random = False
            while n_random <= 10000:
                k = np.random.binomial(n_random, p)
                p_hat_random = k / n_random
                SE_random = np.sqrt(p_hat_random * (1 - p_hat_random) / n_random)
                z_random = (p_hat_random - p0) / SE_random if SE_random > 0 else 0
                if p_hat_random - z_alpha_2 * SE_random > p0:
                    n_random_list.append(n_random)
                    found_random = True
                    break
                n_random += 1
            if not found_random:
                n_random_list.append(10000)  # 未找到
        
        # 计算统计量
        n_avg = np.mean(n_random_list)
        n_std = np.std(n_random_list)
        error_percent = ((n_avg - n_theoretical) / n_theoretical) * 100
        error_range = f"{n_avg-n_std:.0f}-{n_avg+n_std:.0f}"
        
        # 调试信息
        success_rate = (len([x for x in n_random_list if x < 10000]) / len(n_random_list)) * 100
        
        print(f"{p:8.3f} {n_theoretical:8d} {n_avg:10.1f} {n_std:10.1f} {error_percent:8.1f} {error_range:>12} (成功率:{success_rate:.1f}%)")

if __name__ == "__main__":
    # 标称次品率p0=0.1，显著性水平alpha=0.05
    print("=== 原始结果 ===")
    find_min_n_for_p(p0=0.1, alpha=0.05, delta_start=0.001, delta_end=0.01, delta_step=0.001, max_n=10000)
    print("\n")
    calculate_error_range()